/* 
NAME: corruption.do
USING .dta file(s): 
	corruption.dta
	noc.dta
	aid-merge.dta


USING .do file(s): cowcodes.do

DESCRIPTION: This program merges data together, creates variable transformations (including 
	a personalism index), and analyses the relationship between personalism index and FPCA
	corruption cases.

AUTHOR: Joseph Wright
ORIGIN DATE: 11.20.15
LAST UPDATE: 02.16.17

*/

use   "$dir\temp.dta",clear
cd "$dir\corruption"
sort cow year

** merge aid data **
merge cow year using aid-merge
tab _merge
drop if _merge==2
rename _merge _merge2
recode milaid econaid l12milaid lnl12milaid l12econaid lnl12econaid (.=0)
sort cow year
save temp, replace

** merge corruption and personal data **
use corruption, clear
egen tag = tag(country year)
keep if tag==1
drop tag
gen cowcode=.
qui do cowcodes
sort cow year
merge cow year using temp
tab _merge
gen case = _merge==3 if _merge~=1
rename _merge _merge1
gen syear  = year(gwf_start)
egen minyear = min(syear), by(cow)
replace minyear = minyear-1918
replace minyear = 0 if minyear<0 & gwf_case_fail~=.
sort cow  year
save temp, replace

** merge oil company ownership **
use noc, clear
gen cowcode=.
qui do cowcodes
sort cow 
merge cow using temp
tab _merge
rename _merge _merge3
sort cow year
save temp, replace

 
** keep only years after FPCA creation (1977) & oil producing countries **
keep if year>1977 & year<=2010 & _merge3==3   & gwf_duration~=. & oecd1==0
sort cow year
save temp, replace

****************************
** Cross-section analysis **
****************************
	tsset cow year
	egen maxcase = max(case), by(gwf_casename)
	egen meanoil = mean(oilpc),by(gwf_casename)
	egen meanpop = mean(lpop),by(gwf_casename)
	egen meangdp = mean(lgdpcap),by(gwf_casename)
	egen meanmil  =mean(lnl12milaid),by(gwf_casename)
	egen meanecon  =mean(lnl12econaid),by(gwf_casename)
	egen meanpers = mean(pers),by(gwf_casename)

	replace meanecon = meanecon/4
	replace meanmil = meanmil/4
	egen ctag = tag(gwf_casename)
	ttest meanpers if ctag==1, by(maxcase)
	krls maxcase meanpers meanpop meangdp america if ctag==1 
	est store m1
	krls maxcase meanpers meanpop meangdp regnoc meanoil america if ctag==1 
	est store m2
	krls maxcase meanpers meanpop meangdp regnoc meanoil meanecon meanmil america if ctag==1,deriv(m3)
	est store m3
 				global color1="gs1"
				global color2="gs8"
				global color3="gs12"

				global vars = 1 /* number of variables we want to not display */
				gen count=_n
				gen name = ""
				forval m = 1/3 {
					gen hi`m' =.
					gen lo`m' =.
					gen mhi`m'  =.
					gen mlo`m' =.
					gen b`m' =.
					gen v`m'=.   				/* number of variables we want to display */
					qui:est restore m`m' 
					matrix O =  e(Output)
					scalar r = rowsof(O)
					local r =  r
					replace v`m'=r- $vars
					forval c = 1/`r'  {
						local d1 = `c'+4+ 2*`m' 
						local d = `c'
						if `c'> v`m' {
							local d =`d1'
						}
						local rownms: rown O 
						local rowname: word `c' of `rownms'
						replace name = "`rowname'" if count==`d'
						local beta = O[`c',1]
						local var = O[`c',2]
						matrix d==(0,0,0,0,0\0,0,0,0,0)
						matrix d[1,3]=  `beta'
						matrix d[1,5] =  `beta' + 1.96*`var'
						matrix d[1,1] =  `beta' - 1.96*`var'
						matrix d[1,4] =  `beta' + 1.65*`var'
						matrix d[1,2] =  `beta' - 1.65*`var'
						qui replace hi`m'    =  d[1,5] if count==`d'
						qui replace lo`m'    =  d[1,1] if count==`d'
						qui replace mhi`m'    =  d[1,4] if count==`d'
						qui replace mlo`m'    =  d[1,2] if count==`d'
						qui replace b`m'  =  d[1,3] if count==`d'
					}
					sum name hi`m'  lo`m'  b`m' v`m' count
				}
				gen count1=(count*-1)+8
				gen count2=count1-.333
				gen count3=count2-.333
				twoway (scatter count1 b1 if count>=1 & count<=7,ylab(1(1)7,glcolor(gs16)) xlab(0(.1).3) mcolor($color1) msymbol(T) yscale(range(0.75 7.25)) ///
				xtitle(Estimate) xline(0,lpat(dash))) (rspike hi1 lo1 count1 if count>=1 &count<=7, horizontal ytitle("") title(Cross-sectional,size(medium)) ///
				ylab(6.67 "{bf:Personalism}" 5.67 "GDP pc" 4.67 "Population" 3.67 "Oil rents" 2.67 "Regulatory NOC" 1.67 "Economic aid" 0.67 "Military aid") ///
				lcolor($color1) lwidth(medthin) legend(lab(1 "base model") lab(4 "add oil") lab(7 "add aid") order(1 4 7) col(3) pos(6) ring(.5)))		///
				(rspike mhi1 mlo1 count1 if count>=1 &count<=7, lwidth(thick) lcolor($color1) horizontal saving(h1.gph,replace))  ///
				///
				(scatter count2 b2 if count>=1 & count<=7,mcolor($color2) msymbol(S)) ///
				(rspike hi2 lo2 count2 if count>=1 &count<=7, horizontal lcolor($color2) lwidth(medthin)) ///
				(rspike mhi2 mlo2 count2 if count>=1 &count<=7, horizontal lcolor($color2) lwidth(thick)) ///
				///
				(scatter count3 b3 if count>=1 & count<=7,mcolor($color3) msymbol(plus) msize(vlarge)) ///
				(rspike hi3 lo3 count3 if count>=1 &count<=7, horizontal lcolor($color3) lwidth(medthin)) ///
				(rspike mhi3 mlo3 count3 if count>=1 &count<=7, horizontal lcolor($color3) lwidth(thick) ///
				yline(6.165,lpat(solid)lcol(gs15)) yline(5.165,lpat(solid)lcol(gs15)) yline(4.165,lpat(solid)lcol(gs15))  ///
				yline(3.165,lpat(solid)lcol(gs15)) yline(2.165,lpat(solid)lcol(gs15)) yline(1.165,lpat(solid)lcol(gs15)))
				drop name hi* lo* b* count* mhi* mlo* v1 v2 v3
 
**************************
** RE Logits, time-vary **
**************************
	replace lnl12econ = lnl12econ/4
	replace lnl12mil = lnl12mil/4
	gen yr = year
	replace yr = 2005 if yr==2001 | yr==1997 | yr==2003 | yr==2004 | yr==1982 | yr==1983
	xi:xtlogit case i.yr lpop lgdpcap pers   america asia easia ssa meast if yr>1980,vce(cluster cow)
	est store l1
	xi:xtlogit case i.yr lpop lgdpcap pers oilpc regnoc america asia easia ssa meast if yr>1980,vce(cluster cow)
	est store l2
	xi:xtlogit case i.yr lpop lgdpcap pers oilpc regnoc lnl12* america asia easia ssa meast if yr>1980,vce(cluster cow)
	est store l3
	label var lgdpcap "GDP pc (log)"
	label var lpop "Pop (log)"
	label var oilpc "Oil pc (log)"
	label var pers "{bf:Personalism}"
	label var regnoc "Reg NOC"
	label var lnl12mil  "Military aid"
	label var lnl12econ  "Economic aid"
	coefplot (l1, msymbol(T) mcol($color1) ciopts(lpat(solid) lcol($color1 $color1)) ) (l2, msymbol(S) mcol($color2) ciopts(lpat(solid) lcol($color2 $color2)) ) /*
		*/ (l3, msymbol(plus)mcol($color3) ciopts(lpat(solid) lcol($color3 $color3)) ), title("RE logit", size(medsmall) height(5))/*
		*/ scheme(lean2) drop(_cons _Iy* americas asia easia ssa meast) order (pers lpop lgdpcap oilpc regnoc lnl12econ lnl12mil) /*
		*/ xlab(0 (2) 8) xline(0, lpattern(dash)) grid(glcolor(gs15)) mfcolor(white) /*
		*/ ysize(2) xsize(2.5)   /*
		*/ legend(label(3 "base model") label(6 "add oil") label(9 "add aid") size(small) pos(6) ring(1.5) col(3))  /*
		*/ levels(95 90) xtitle("  Coefficient estimate", height(6)size(small))
	graph export "$dir\golden\Corruption-RElogit.pdf", as(pdf)   replace
	 xtlogit case period* lpop lgdp pers oilpc regnoc america asia easia ssa meast,vce(cluster cow)
 
************************
** Conditional logits **
************************
	gen time = year-1978
	gen xlpop =lpop*10 
	clogit case time xlpop lgdpcap pers if year>1980,group(cow)  
	est store cl1
	clogit case time xlpop lgdpcap pers oilpc if year>1980,group(cow)
	est store cl2
	clogit case time xlpop lgdpcap pers oilpc lnl12* if year>1980,group(cow)   
	est store cl3
	label var time "Time trend"
	label var xlpop "Population"
	coefplot (cl1, msymbol(T) mcol($color1) ciopts(lpat(solid) lcol($color1 $color1)) ) (cl2, msymbol(S) mcol($color2) ciopts(lpat(solid) lcol($color2 $color2)) ) /*
		*/ (cl3, msymbol(plus)mcol($color3) ciopts(lpat(solid) lcol($color3 $color3)) ), title("Conditional logit", size(medsmall) height(5))/*
		*/ scheme(lean2) drop(_cons _Iy* americas asia easia ssa meast) order (pers time xlpop lgdpcap oilpc lnl12econ lnl12mil) /*
		*/ xlab(-4 (2) 6) xline(0, lpattern(dash)) grid(glcolor(gs15)) mfcolor(white) /*
		*/ ysize(2) xsize(2.5)   /*
		*/ legend(label(3 "base model") label(6 "add oil") label(9 "add aid") size(small) pos(6) ring(1.5) col(3))  /*
		*/ levels(95 90) xtitle("  Coefficient estimate", height(6)size(small))
	graph export "$dir\golden\Corruption-Clogit.pdf", as(pdf)   replace
	
	
*****************
** Firth logit **
*****************
* post 1980 *
firthlogit case regnoc lpop lgdpcap pers oilpc america asia easia ssa meast if year>1978
* add reg NOC *
 firthlogit case regnoc lpop lgdpcap pers oilpc america asia easia ssa meast if  year>1978

********************
** KRLS time-vary **
********************
	 krls case pers lgdpcap lpop time america ssa meast asia easia if year>1996,deriv(k1) vcov
	est store c1
	 krls case pers lgdpcap lpop oilpc regnoc time america ssa meast asia easia if year>1996,deriv(k2) vcov
	est store c2
	 krls case pers lgdpcap lpop oilpc regnoc lnl12econ lnl12mil time america ssa meast asia easia if year>1996,deriv(k3) vcov
	 est store c3
	* Check with missing lgdpcap data filled-in for Argentina *
	qui:reg lgdpcap lgdp lpop i.cow i.year
	predict hat 
	 krls case pers hat lpop oilpc regnoc lnl12econ lnl12mil time america ssa meast asia easia if year>1996,deriv(k4) vcov
		twoway (kdensity k3_pers  , lcolor($color1) xline(0,lpattern(dash)) ///
		bw(0.015) xtitle(Derivatives) ylab(0 (2) 6,glcolor(gs15)))  (kdensity k4_pers,  ///
		lcolor($color2) bw(.015) legend(label(1 "n=615") lab(2 "n=619") lab(3 "",,off) lab(4 "",off) ///
		order(1 2)  pos(2) ring(0) col(1)) ytitle(Density) )
		
	 * Other time frames *
	 krls case pers lgdpcap lpop oilpc regnoc time america ssa meast asia easia if year>1978, vcov
	 krls case pers lgdpcap lpop oilpc regnoc time america ssa meast asia easia if year>2000, vcov
	 krls case pers lgdpcap lpop oilpc regnoc      america ssa meast asia easia if year>2006, vcov			
	
	* KRLS estimates *
		global color1="gs1"
		global color2="gs8"
		global color3="gs12"
				global vars = 6 /* number of variables we want to not display */
				gen count=_n
				gen name = ""
				forval m = 1/3 {
					gen hi`m' =.
					gen lo`m' =.
					gen mhi`m'  =.
					gen mlo`m' =.
					gen b`m' =.
					gen v`m'=.   				/* number of variables we want to display */
					qui:est restore c`m' 
					matrix O =  e(Output)
					scalar r = rowsof(O)
					local r =  r
					replace v`m'=r- $vars
					forval c = 1/`r'  {
						local d1 = `c'+2+ 2*`m' 
						local d = `c'
						if `c'> v`m' {
							local d =`d1'
						}
						local rownms: rown O 
						local rowname: word `c' of `rownms'
						replace name = "`rowname'" if count==`d'
						local beta = O[`c',1]
						local var = O[`c',2]
						matrix d==(0,0,0,0,0\0,0,0,0,0)
						matrix d[1,3]=  `beta'
						matrix d[1,5] =  `beta' + 1.96*`var'
						matrix d[1,1] =  `beta' - 1.96*`var'
						matrix d[1,4] =  `beta' + 1.65*`var'
						matrix d[1,2] =  `beta' - 1.65*`var'
						qui replace hi`m'    =  d[1,5] if count==`d'
						qui replace lo`m'    =  d[1,1] if count==`d'
						qui replace mhi`m'    =  d[1,4] if count==`d'
						qui replace mlo`m'    =  d[1,2] if count==`d'
						qui replace b`m'  =  d[1,3] if count==`d'
					}
					sum name hi`m'  lo`m'  b`m' v`m' count
				}
				gen count1=(count*-1)+8
				gen count2=count1-.333
				gen count3=count2-.333
	 
				twoway (scatter count1 b1 if count>=1 & count<=7,ylab(1(1)7,glcolor(gs16)) mcolor($color1) msymbol(T) yscale(range(0.75 7.25)) ///
				xtitle(Estimate,size(medsmall)) xline(0,lpat(dash))) (rspike hi1 lo1 count1 if count>=1 &count<=7, horizontal ytitle("") title(Time-varying,size(medium)) ///
				ylab(6.67 "{bf:Personalism}" 5.67 "GDP pc" 4.67 "Population" 3.67 "Oil rents" 2.67 "Regulatory NOC" 1.67 "Economic aid" 0.67 "Military aid") ///
				lcolor($color1) lwidth(medthin) legend(lab(1 "base model") lab(4 "add oil") lab(7 "add aid") order(1 4 7) col(3) pos(6) ring(.5)))		///
				(rspike mhi1 mlo1 count1 if count>=1 &count<=7, lwidth(thick) lcolor($color1) horizontal saving(h2.gph,replace))  ///
				///
				(scatter count2 b2 if count>=1 & count<=7,mcolor($color2) msymbol(S)) ///
				(rspike hi2 lo2 count2 if count>=1 &count<=7, horizontal lcolor($color2) lwidth(medthin)) ///
				(rspike mhi2 mlo2 count2 if count>=1 &count<=7, horizontal lcolor($color2) lwidth(thick)) ///
				///
				(scatter count3 b3 if count>=1 & count<=7,mcolor($color3) msymbol(plus) msize(vlarge)) ///
				(rspike hi3 lo3 count3 if count>=1 &count<=7, horizontal lcolor($color3) lwidth(medthin)) ///
				(rspike mhi3 mlo3 count3 if count>=1 &count<=7, horizontal lcolor($color3) lwidth(thick) ///
				yline(6.165,lpat(solid)lcol(gs15)) yline(5.165,lpat(solid)lcol(gs15)) yline(4.165,lpat(solid)lcol(gs15))  ///
				yline(3.165,lpat(solid)lcol(gs15)) yline(2.165,lpat(solid)lcol(gs15)) yline(1.165,lpat(solid)lcol(gs15)))
				drop name hi* lo* b* count* mhi* mlo* v1 v2 v3
				gr combine h1.gph h2.gph,xsize(4) ysize(2.25)
				graph export "$dir\golden\Corruption-KRLS.pdf", as(pdf)   replace
				graph export "$dir\golden\ISQ-Figure-10.png", as(png)   replace
 
		* KRLS densities, by reg NOC *
		twoway (kdensity k3_pers if regnoc==1, lcolor($color1) xline(0,lpattern(dash)) ///
		bw(0.015) xtitle(Derivatives) ylab(0 (2) 6,glcolor(gs15)))  (kdensity k3_pers if regnoc==0,  ///
		lcolor($color2) bw(.015) legend(label(1 "Reg NOC") lab(2 "no Reg NOC") lab(3 "",,off) lab(4 "",off) ///
		order(1 2)  pos(2) ring(0) col(1)) ytitle(Density)  ///
		text(2.3  -.185 " {bf:{&mu} =0.011}",linegap(-1.3)place(n)) ///
		text(5.25  .15 " {bf:{&mu} =0.093}",linegap(-1.3)place(n)))  ///
		(pcarrowi 2.25  -.175 1.6  -.09)  ///
		(pcarrowi 5.25  .15 4.5  .07,saving(h2,replace) title(Model 3))

		ttest k2_pers,by(regnoc)
		twoway (kdensity k2_pers if regnoc==1, lcolor($color1) xline(0,lpattern(dash)) ///
		bw(0.015) xtitle(Derivatives) ylab(0 (2) 6,glcolor(gs15)))  (kdensity k2_pers if regnoc==0,  ///
		lcolor($color2) bw(.015) legend(label(1 "Reg NOC") lab(2 "no Reg NOC") lab(3 "",,off) lab(4 "",off) ///
		order(1 2)  pos(2) ring(0) col(1)) ytitle(Density)  ///
		text(1.8  .38 " {bf:{&mu} =0.132}",linegap(-1.3)place(n)) ///
		text(5.5  .15 " {bf:{&mu} =0.038}",linegap(-1.3)place(n)))  ///
		(pcarrowi 1.75  .38 1.31  .335)  ///
		(pcarrowi 5.5  .15 5.2  .07,saving(h1,replace) title(Model 2))
		gr combine h1.gph h2.gph,xsize(4.6) ysize(2.25)
		graph export "$dir\golden\Corruption-KRLS-Densities.pdf", as(pdf)   replace
		
  erase h1.gph 
  erase h2.gph

  ****************************************** End of Corruption file *************************************
